
%%% The Ambiguity of Fishing for Fun, October 2022
%%% This file calibrates, for each of the four utility models considered, the total number of choice occasions so that the number of expected
%%% trips and landings match the actual data for the state of New Jersey in 2019

%Loading input files
%Cost per trip and demographic information
results=importdata('ReadyDataFluke.txt'); 
cost=importdata('trip_costs_NE.xlsx'); 
cost=cost.data(cost.textdata(:,2)=="NJ");
C=unique(cost)/1000; 
[r5,c5]=find(strcmp(results.textdata,'Age'));
[r6,c6]=find(strcmp(results.textdata,'Male'));
y=results.data(:,c5:c6);
Y=unique(y,'rows');
clearvars results
%Catch-per-trip
CatchPerTrip=importdata('CatchPerTrip2019.mat');
SFreshaffle=CatchPerTrip.SFreshaffle2019;
Sreshaffle=CatchPerTrip.Sreshaffle2019;
BSBreshaffle=CatchPerTrip.BSBreshaffle2019;
Wreshaffle=CatchPerTrip.Wreshaffle2019;
RDreshaffle=CatchPerTrip.RDreshaffle2019;
clearvars SFreshaffle19 Sreshaffle19 BSBreshaffle19 Wreshaffle19 RDreshaffle19
%Catch-at-length
PSFSizes=importdata('ProbFlukeSizes2019.txt'); 
SFpossibleSizes=PSFSizes(:,1);
PSSizes=importdata('ProbScupSizes2019.txt');
SpossibleSizes=PSSizes(:,1);
PBSBSizes=importdata('ProbBSBSizes2019.txt'); 
BSBpossibleSizes=PBSBSizes(:,1);
catchSF_max=CatchPerTrip.catchSF_max19;
catchSF_min=CatchPerTrip.catchSF_min19;
%Parameter estimates for the four utility models
parameters_cara=importdata('parameters_cara.txt');
parameters_ambiguity=importdata('parameters_ambiguity.txt'); 
parameters_mixed=importdata('parameters_mixed.xlsx'); 
paramLinear=importdata('paramlinear.txt'); 
sigma_cara=importdata('sigma_cara.txt');
sigma_ambiguity=importdata('sigma_ambiguity.txt'); 
sigma_mixed=importdata('sigma_mixed.xlsx'); 
sigmaLinear=importdata('sigmaL.txt');

%Regulations
SFbaglimit=3;  SFsizelimit=18; Sbaglimit=50; Ssizelimit=9; BSBbaglimit=10; BSBsizelimit=12.5; Wbaglimit=1; Wsizelimit=13;

%Actual number of trips in New Jersey in 2019
NETrips=importdata('G:\My Drive\RESEARCH\WORKING PAPERS\McConnell_Holzer on Responses to Risky Catch\Matlab code\simulations\data from Lou\mrip_directed_sf_trips_2009_2020.xlsx'); 
NETrips.textdata=NETrips.textdata(26:end,:);
[r7,c7]=find(contains(NETrips.textdata,'NEW JERSEY'));
NETrips.data=horzcat(NETrips.data(:,1),NETrips.data(:,3));
NETrips.data=NETrips.data(r7,:);
[r8,c8]=find(NETrips.data==2019);
NETrips.data=NETrips.data(r8,2);
TRIPS=round(NETrips.data/1000);
clearvars NETrips

%Observed average keep per trip for the differernt species
CUTOFFS=importdata('AvgeCatchTrip.xlsx'); 
[r7,c7]=find(contains(CUTOFFS.textdata,'NJ'));
[r8,c8]=find(contains(CUTOFFS.textdata,'summerflounder'));
[r9,c9]=find(contains(CUTOFFS.textdata,'blackseabass'));
[r10,c10]=find(contains(CUTOFFS.textdata,'scup'));
[r11,c11]=find(contains(CUTOFFS.textdata,'weakfish'));
[r12,c12]=find(contains(CUTOFFS.textdata,'reddrum'));
[r13,c13]=find(CUTOFFS.data==2019);

cutoff_sf=CUTOFFS.data(intersect(intersect(r7,r8),r13),6); 
cutoff_bsb=CUTOFFS.data(intersect(intersect(r7,r9),r13),6);
cutoff_s=CUTOFFS.data(intersect(intersect(r7,r10),r13),6);
cutoff_w=CUTOFFS.data(intersect(intersect(r7,r11),r13),6);
cutoff_rd=CUTOFFS.data(intersect(intersect(r7,r12),r13),6);

T=100; %Setting number of parameters' draws to construct confidence intervals
M=25000;

%Initializing vectors
[AggregateKeepSF_mixed,AggregateKeepSF_cara,AggregateKeepSF_ambiguity,AggregateKeepSFLinear,AggregateReleaseSF_mixed, AggregateReleaseSF_cara,AggregateReleaseSF_ambiguity,AggregateReleaseSFLinear,AggregateKeepSnumber_cara,AggregateKeepBSBnumber_cara,AggregateReleasedBSBnumber_cara,AggregateKeepWLinear,AggregateKeepWnumberLinear,AggregateReleaseWLinear,AggregateReleasedWnumberLinear,AggregateKeepBSB_cara,AggregateReleaseSLinear,AggregateReleaseBSB_mixed,AggregateKeepSFnumber_mixed]=deal(zeros(M,T));
[AggregateKeepS_mixed,AggregateKeepS_cara, AggregateKeepS_ambiguity,AggregateKeepSLinear,AggregateReleaseS_mixed,AggregateReleaseS_cara,AggregateReleaseS_ambiguity,AggregateReleaseSLinearAggregateKeepBSB_cara,AggregateKeepBSB_mixed,AggregateKeepBSB_ambiguity,AggregateKeepBSBLinear,AggregateReleaseBSB_cara,AggregateReleaseBSB,AggregateReleaseBSB_ambiguity,AggregateReleaseBSBLinear,AggregateKeepSFnumber_cara,AggregateReleasedSnumber_cara,AggregateReleasedSFnumber_cara,AggregateKeepSFnumber_mixe,AggregateKeepSFnumber_ambiguity,AggregateKeepSnumber_ambiguity,AggregateReleasedSnumber_mixed]=deal(zeros(M,T));
[AggregateReleasedSFnumber_mixed,AggregateReleaseSFnumber_ambiguity,AggregateKeepSFnumberLinear,AggregateReleasedSFnumberLinear,AggregateKeepSnumber_mixed,AggregateReleasedSnumber_ambiguity,AggregateKeepSnumberLinear,AggregateReleasedBSBnumberLinear,AggregateReleasedBSBnumber_mixed,AggregateKeepBSBnumberLinear,AggregateReleasedSnumberLinear,AggregateKeepBSBnumber_ambiguity,AggregateReleasedBSBnumber_ambiguity,AggregateKeepBSBnumber_mixed,TotalWelfare_mixed,TotalWelfare_cara,TotalWelfare_ambiguity,TotalWelfareLinear,NumberOftrips_cara,NumberOftrips_mixed,NumberOftrips_ambiguity,NumberOftripsLinear]=deal(zeros(M,T));
[AggregateKeepW_ambiguity,AggregateKeepWnumber_ambiguity,AggregateReleaseW_ambiguity,AggregateReleasedWnumber_ambiguity,AggregateKeepW_mixed,AggregateKeepWnumber_mixed,AggregateReleaseW_mixed,AggregateReleasedWnumber_mixed,AggregateKeepW_cara,AggregateKeepWnumber_cara,AggregateReleaseW_cara,AggregateReleasedWnumber_cara]=deal(zeros(M,T));
[vartheta_cara,vartheta_ambiguity,varthetaLinear]=deal(zeros(T,7)); 

tic; 
rng(4560); 
N=30;  %Number of draws for each potential trip

%initializing vectors
[u_mixed,u_cara,uLinear,u_ambiguity]=deal(zeros(T,N));
[v_mixed,v_cara,vLinear,v_ambiguity,opt_out_mixed,opt_out_ambiguity,opt_out_cara,probTripLinear,opt_outLinear,probTripCARA,probTripAmbiguity,probTripMixed]=deal(zeros(M,T));
[SFkept,SFreleased,SFkeptLbs,Skept,Sreleased,SkeptLbs,Wkept,Wreleased,WkeptLbs,WreleasedLbs,SFreleasedLbs,SreleasedLbs,BSBreleasedLbs,BSBkept,BSBreleased,BSBkeptLbs,SFkeptNumbers]=deal(zeros(N,T));
[Wnumber,ReleasedSFnumber,ReleasedSnumber,ReleasedBSBnumber,BSBnumber,Snumber,SFnumber,SFLbs,SLbs,BSBLbs,WLb,ReleasedWnumber,ReleasedBSBLbs,ReleasedWLbs,WLbs,ReleasedSFLbs,ReleasedSLbs]=deal(zeros(M,T)); 
par_mixed=zeros(T,(length(parameters_mixed)-1)/2+1,M);

myObj=@(p)pstar(M,N,SFreshaffle(:,:,1),p,cutoff_sf,SFbaglimit);
options=optimset('MaxFunEvals',50000,'TolX',1e-15,'MaxIter',10000, 'Display','on');
[p_cutoff_SF,fval,exitflag,output]=fzero(myObj,0.9,options);

%Probability the size limit binds for each species given MRIP data
myObj=@(p)pstar(M,N,BSBreshaffle(:,:,1),p,cutoff_bsb,BSBbaglimit);
[p_cutoff_BSB,fval_BSB,exitflag_BSBS,output_BSB]=fzero(myObj,0.5,options);

if cutoff_s~=0
myObj=@(p)pstar(M,N,Sreshaffle(:,:,1),p,cutoff_s,Sbaglimit);
[p_cutoff_S,fval_S,exitflag_S,output_S]=fzero(myObj,0.2,options);
end

if cutoff_w>0.05
myObj=@(p)pstar(M,N,Wreshaffle(:,:,1),p,cutoff_w,Wbaglimit);
[p_cutoff_W,fval_W,exitflag_W,output_W]=fzero(myObj,0.9,options);
else
p_cutoff_W=1;    
end

%Drawing parameters from estimates
parLinear=mvnrnd(paramLinear,sigmaLinear,T);
par_cara=mvnrnd(parameters_cara,sigma_cara,T);
par_ambiguity=mvnrnd(parameters_ambiguity,sigma_ambiguity,T);
param_mixed=mvnrnd(parameters_mixed,sigma_mixed,T);

%1. Alpha-maxmin Utility Model
  k=0;
  while (k<M) && (mean(sum(NumberOftrips_ambiguity,1))<=TRIPS)     
    k=k+1;
    %Random draw for trip non-catch characteristics
    IndvCharac=Y(randperm(end),:); 
    CostCharac=datasample(C,1); 
    CostCharacOpt_out=datasample(C,1);
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    SEncounters=Sreshaffle(k,:,t); 
    BSBEncounters=BSBreshaffle(k,:,t);       
    WEncounters=Wreshaffle(k,:,t);   
    RDEncounters=RDreshaffle(k,:,t);  
    
    for j=1:N   
     %Random draws for summer flounder
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
             
     %Random draws for scup
     if SEncounters(j)>0
     z=rand(1,SEncounters(j));
     Slegal=(z>=p_cutoff_S);
     Skept(j,t)=min(Sbaglimit,sum(Slegal,'all')); 
     Sreleased(j,t)=SEncounters(j)-Skept(j,t);       
     else
     Skept(j,t)=0;
     Sreleased(j,t)=0;
     end
     
     %Random draws for black sea bass
     if BSBEncounters(j)>0
     z=rand(1,BSBEncounters(j));
     BSBlegal=(z>=p_cutoff_BSB);
     BSBkept(j,t)=min(BSBbaglimit,sum(BSBlegal,'all'));
     BSBreleased(j,t)=BSBEncounters(j)-BSBkept(j,t);       
     else
     BSBkept(j,t)=0;
     BSBreleased(j,t)=0;
     end
  
     %Random draws for weakfish
     if WEncounters(j)>0
     z=rand(1,WEncounters(j));
     Wlegal=(z>=p_cutoff_W);
     Wkept(j,t)=min(Wbaglimit,sum(Wlegal,'all'));
     Wreleased(j,t)=WEncounters(j)-Wkept(j,t);       
     else
     Wkept(j,t)=0;
     Wreleased(j,t)=0;
     end
     
     u_ambiguity(t,j)=par_ambiguity(t,15)*(1-exp(-par_ambiguity(t,5)*(par_ambiguity(t,2)*0/100+par_ambiguity(t,1)*0/100+par_ambiguity(t,4)*(Skept(j,t)+BSBkept(j,t)+Wkept(j,t))/100+par_ambiguity(t,3)*(Sreleased(j,t)+BSBreleased(j,t)+Wreleased(j,t))/100)))/par_ambiguity(t,5)+(1-par_ambiguity(t,15))*(1-exp(-par_ambiguity(t,5)*(par_ambiguity(t,2)*SFbaglimit/100+par_ambiguity(t,1)*(catchSF_max-SFbaglimit)/100+par_ambiguity(t,4)*(Skept(j,t)+BSBkept(j,t)+Wkept(j,t))/100+par_ambiguity(t,3)*(Sreleased(j,t)+BSBreleased(j,t)+Wreleased(j,t))/100)))/par_ambiguity(t,5);   
     end %End of the N random draws per individual
    
    vartheta_ambiguity(t,:)=[par_ambiguity(t,8),par_ambiguity(t,9),par_ambiguity(t,10),par_ambiguity(t,11),par_ambiguity(t,12),par_ambiguity(t,13), par_ambiguity(t,14)];
    
    %Calculating the individual utility for individual k of M
    v_ambiguity(k,t)=mean(u_ambiguity(t,:))+par_ambiguity(t,6)*CostCharac(1);
    opt_out_ambiguity(k,t)=par_ambiguity(t,7)+IndvCharac(1,:)*vartheta_ambiguity(t,:)'+par_ambiguity(t,6)*CostCharacOpt_out(1);
    
    %Calculating the probabilities of taking the trip
    probTripAmbiguity(k,t)=exp(v_ambiguity(k,t))/(exp(v_ambiguity(k,t))+exp(opt_out_ambiguity(k,t)));
    
    %Calculating catch per trip
    SFLbs(k,t)=mean(SFkeptLbs(:,t));
    SFnumber(k,t)=mean(SFkept(:,t));
    SLbs(k,t)=mean(SkeptLbs(:,t));
    Snumber(k,t)=mean(Skept(:,t));
    BSBLbs(k,t)=mean(BSBkeptLbs(:,t));
    BSBnumber(k,t)=mean(BSBkept(:,t));
    WLbs(k,t)=mean(WkeptLbs(:,t));
    Wnumber(k,t)=mean(Wkept(:,t));
    ReleasedSFLbs(k,t)=mean(SFreleasedLbs(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    ReleasedSLbs(k,t)=mean(SreleasedLbs(:,t));
    ReleasedSnumber(k,t)=mean(Sreleased(:,t));
    ReleasedBSBLbs(k,t)=mean(BSBreleasedLbs(:,t));
    ReleasedBSBnumber(k,t)=mean(BSBreleased(:,t));
    ReleasedWLbs(k,t)=mean(WreleasedLbs(:,t));
    ReleasedWnumber(k,t)=mean(Wreleased(:,t));
    
    %Calculating total number of trips, aggregate catch and welfare
    AggregateKeepSF_ambiguity(k,t)=probTripAmbiguity(k,t)*SFLbs(k,t);
    AggregateKeepSFnumber_ambiguity(k,t)=probTripAmbiguity(k,t)*SFnumber(k,t);
    AggregateReleaseSF_ambiguity(k,t)=probTripAmbiguity(k,t)*ReleasedSFLbs(k,t);
    AggregateReleaseSFnumber_ambiguity(k,t)=probTripAmbiguity(k,t)*ReleasedSFnumber(k,t);
    
    AggregateKeepS_ambiguity(k,t)=probTripAmbiguity(k,t)*SLbs(k,t);
    AggregateKeepSnumber_ambiguity(k,t)=probTripAmbiguity(k,t)*Snumber(k,t);
    AggregateReleaseS_ambiguity(k,t)=probTripAmbiguity(k,t)*ReleasedSLbs(k,t);
    AggregateReleasedSnumber_ambiguity(k,t)=probTripAmbiguity(k,t)*ReleasedSnumber(k,t);
    
    AggregateKeepBSB_ambiguity(k,t)=probTripAmbiguity(k,t)*BSBLbs(k,t);
    AggregateKeepBSBnumber_ambiguity(k,t)=probTripAmbiguity(k,t)*BSBnumber(k,t);
    AggregateReleaseBSB_ambiguity(k,t)=probTripAmbiguity(k,t)*ReleasedBSBLbs(k,t); 
    AggregateReleasedBSBnumber_ambiguity(k,t)=probTripAmbiguity(k,t)*ReleasedBSBnumber(k,t);
    
    AggregateKeepW_ambiguity(k,t)=probTripAmbiguity(k,t)*WLbs(k,t);
    AggregateKeepWnumber_ambiguity(k,t)=probTripAmbiguity(k,t)*Wnumber(k,t);
    AggregateReleaseW_ambiguity(k,t)=probTripAmbiguity(k,t)*ReleasedWLbs(k,t); 
    AggregateReleasedWnumber_ambiguity(k,t)=probTripAmbiguity(k,t)*ReleasedWnumber(k,t);
        
    TotalWelfare_ambiguity(k,t)=log(exp(v_ambiguity(k,t))+exp(opt_out_ambiguity(k,t)))/(-par_ambiguity(t,6)/1000);
    NumberOftrips_ambiguity(k,t)=probTripAmbiguity(k,t);
    end
  end

k1=k;
  
  %2. Linear Utility Model (mlogit estimator)
    k=0;
   while (k<M) && (mean(sum(NumberOftrips_mixed,1))<=TRIPS) 
    k=k+1;
    %Random draw for trip non-catch characteristics
    IndvCharac=Y(randperm(end),:); 
    CostCharac=datasample(C,1); 
    CostCharacOpt_out=datasample(C,1);
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    SEncounters=Sreshaffle(k,:,t); 
    BSBEncounters=BSBreshaffle(k,:,t);    
    WEncounters=Wreshaffle(k,:,t); 

    for j=1:N   
     %Random draws for summer flounder
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
             
     %Random draws for scup
     if SEncounters(j)>0
     z=rand(1,SEncounters(j));
     Slegal=(z>=p_cutoff_S);
     Skept(j,t)=sum(Slegal,'all'); 
     Sreleased(j,t)=SEncounters(j)-Skept(j,t);       
     else
     Skept(j,t)=0;
     Sreleased(j,t)=0;
     end
     
     %Random draws for black sea bass
     if BSBEncounters(j)>0
     z=rand(1,BSBEncounters(j));
     BSBlegal=(z>=p_cutoff_BSB);
     BSBkept(j,t)=min(BSBbaglimit,sum(BSBlegal,'all')); .
     BSBreleased(j,t)=BSBEncounters(j)-BSBkept(j,t);       
     else
     BSBkept(j,t)=0;
     BSBreleased(j,t)=0;
     end
        
     %Random draws for weakfish
     if WEncounters(j)>0
     z=rand(1,WEncounters(j));
     Wlegal=(z>=p_cutoff_W);
     Wkept(j,t)=min(Wbaglimit,sum(Wlegal,'all')); 
     Wreleased(j,t)=WEncounters(j)-Wkept(j,t);       
     else
     Wkept(j,t)=0;
     Wreleased(j,t)=0;
     end   
     
     par_mixed(t,:,k)=[param_mixed(t,1);normrnd(param_mixed(t,2),abs(param_mixed(t,7)));normrnd(param_mixed(t,3),abs(param_mixed(t,8)));normrnd(param_mixed(t,4),abs(param_mixed(t,9)));normrnd(param_mixed(t,5),abs(param_mixed(t,10)));normrnd(param_mixed(t,6),abs(param_mixed(t,11)))]; 
     u_mixed(t,j)=par_mixed(t,3,k)*SFkept(j,t)/100+par_mixed(t,2,k)*SFreleased(j,t)/100+par_mixed(t,5,k)*(Skept(j,t)+BSBkept(j,t)+Wkept(j,t))/100+par_mixed(t,4,k)*(Sreleased(j,t)+BSBreleased(j,t)+Wreleased(j,t))/100;     
     end %End of the N random draws per individual
    
    v_mixed(k,t)=mean(u_mixed(t,:))+par_mixed(t,1,k)*CostCharac(1);
    opt_out_mixed(k,t)=par_mixed(t,6,k);
    
    %Calculating the probabilities of taking the trip
    probTripMixed(k,t)=exp(v_mixed(k,t))/(exp(v_mixed(k,t))+exp(opt_out_mixed(k,t)));
    
    %Calculating catch per trip
    SFLbs(k,t)=mean(SFkeptLbs(:,t));
    SFnumber(k,t)=mean(SFkept(:,t));
    SLbs(k,t)=mean(SkeptLbs(:,t));
    Snumber(k,t)=mean(Skept(:,t));
    BSBLbs(k,t)=mean(BSBkeptLbs(:,t));
    BSBnumber(k,t)=mean(BSBkept(:,t));
    WLbs(k,t)=mean(WkeptLbs(:,t));
    Wnumber(k,t)=mean(Wkept(:,t));
    ReleasedSFLbs(k,t)=mean(SFreleasedLbs(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    ReleasedSLbs(k,t)=mean(SreleasedLbs(:,t));
    ReleasedSnumber(k,t)=mean(Sreleased(:,t));
    ReleasedBSBLbs(k,t)=mean(BSBreleasedLbs(:,t));
    ReleasedBSBnumber(k,t)=mean(BSBreleased(:,t));
    ReleasedWLbs(k,t)=mean(WreleasedLbs(:,t));
    
    %Calculating total number of trips, aggregate catch and welfare:
    AggregateKeepSF_mixed(k,t)=probTripMixed(k,t)*SFLbs(k,t);
    AggregateKeepSFnumber_mixed(k,t)=probTripMixed(k,t)*SFnumber(k,t);
    AggregateReleaseSF_mixed(k,t)=probTripMixed(k,t)*ReleasedSFLbs(k,t);
    AggregateReleasedSFnumber_mixed(k,t)=probTripMixed(k,t)*ReleasedSFnumber(k,t);
    
    AggregateKeepS_mixed(k,t)=probTripMixed(k,t)*SLbs(k,t);
    AggregateKeepSnumber_mixed(k,t)=probTripMixed(k,t)*Snumber(k,t);
    AggregateReleaseS_mixed(k,t)=probTripMixed(k,t)*ReleasedSLbs(k,t);
    AggregateReleasedSnumber_mixed(k,t)=probTripMixed(k,t)*ReleasedSnumber(k,t);
    
    AggregateKeepBSB_mixed(k,t)=probTripMixed(k,t)*BSBLbs(k,t);
    AggregateKeepBSBnumber_mixed(k,t)=probTripMixed(k,t)*BSBnumber(k,t);
    AggregateReleaseBSB_mixed(k,t)=probTripMixed(k,t)*ReleasedBSBLbs(k,t);
    AggregateReleasedBSBnumber_mixed(k,t)=probTripMixed(k,t)*ReleasedBSBnumber(k,t);
    
    AggregateKeepW_mixed(k,t)=probTripMixed(k,t)*WLbs(k,t);
    AggregateKeepWnumber_mixed(k,t)=probTripMixed(k,t)*Wnumber(k,t);
    AggregateReleaseW_mixed(k,t)=probTripMixed(k,t)*ReleasedWLbs(k,t);
    AggregateReleasedWnumber_mixed(k,t)=probTripMixed(k,t)*ReleasedWnumber(k,t);
        
    TotalWelfare_mixed(k,t)=log(exp(v_mixed(k,t))+exp(opt_out_mixed(k,t)))/(-par_mixed(t,1,k)/1000);    
    NumberOftrips_mixed(k,t)=probTripMixed(k,t);
    end
  end

k2=k;

%3. Linear Utility Model (clogit estimator)
    k=0;
  while (k<M) && (mean(sum(NumberOftripsLinear,1))<=TRIPS)  
    k=k+1;
    %Random draw for trip non-catch characteristics
    IndvCharac=Y(randperm(end),:); 
    CostCharac=datasample(C,1); 
    CostCharacOpt_out=datasample(C,1);
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    SEncounters=Sreshaffle(k,:,t); 
    BSBEncounters=BSBreshaffle(k,:,t);      
    WEncounters=Wreshaffle(k,:,t); 

    for j=1:N   
     %Random draws for summer flounder
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
             
     %Random draws for scup
     if SEncounters(j)>0
     z=rand(1,SEncounters(j));
     Slegal=(z>=p_cutoff_S);
     Skept(j,t)=sum(Slegal,'all');
     Sreleased(j,t)=SEncounters(j)-Skept(j,t);       
     else
     Skept(j,t)=0;
     Sreleased(j,t)=0;
     end
     
     %Random draws for black sea bass
     if BSBEncounters(j)>0
     z=rand(1,BSBEncounters(j));
     BSBlegal=(z>=p_cutoff_BSB);
     BSBkept(j,t)=min(BSBbaglimit,sum(BSBlegal,'all'));
     BSBreleased(j,t)=BSBEncounters(j)-BSBkept(j,t);       
     else
     BSBkept(j,t)=0;
     BSBreleased(j,t)=0;
     end
     
     %Random draws for weakfish
     if WEncounters(j)>0
     z=rand(1,WEncounters(j));
     Wlegal=(z>=p_cutoff_W);
     Wkept(j,t)=min(Wbaglimit,sum(Wlegal,'all')); 
     Wreleased(j,t)=WEncounters(j)-Wkept(j,t);       
     else
     Wkept(j,t)=0;
     Wreleased(j,t)=0;
     end
     
     uLinear(t,j)=parLinear(t,2)*SFkept(j,t)/100+parLinear(t,1)*SFreleased(j,t)/100+parLinear(t,4)*(Skept(j,t)+BSBkept(j,t)+Wkept(j,t))/100+parLinear(t,3)*(Sreleased(j,t)+BSBreleased(j,t)+Wreleased(j,t))/100;
     end %End of the N random draws per individual
    
    varthetaLinear(t,:)=[parLinear(t,7),parLinear(t,8), parLinear(t,9), parLinear(t,10), parLinear(t,11),parLinear(t,12), parLinear(t,13)]; 
    
    %Calculating the individual utility for individual k of M
    vLinear(k,t)=mean(uLinear(t,:))+parLinear(t,5)*CostCharac(1); 
    opt_outLinear(k,t)=parLinear(t,6)+IndvCharac(1,:)*varthetaLinear(t,:)'+parLinear(t,5)*CostCharacOpt_out(1);
 
    %Calculating the probabilities of taking the trip
     probTripLinear(k,t)=exp(vLinear(k,t))/(exp(vLinear(k,t))+exp(opt_outLinear(k,t)));
     
    %Calculating catch per trip
    SFLbs(k,t)=mean(SFkeptLbs(:,t));
    SFnumber(k,t)=mean(SFkept(:,t));
    SLbs(k,t)=mean(SkeptLbs(:,t));
    Snumber(k,t)=mean(Skept(:,t));
    BSBLbs(k,t)=mean(BSBkeptLbs(:,t));
    BSBnumber(k,t)=mean(BSBkept(:,t));
    WLbs(k,t)=mean(WkeptLbs(:,t));
    Wnumber(k,t)=mean(Wkept(:,t));
    ReleasedSFLbs(k,t)=mean(SFreleasedLbs(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    ReleasedSLbs(k,t)=mean(SreleasedLbs(:,t));
    ReleasedSnumber(k,t)=mean(Sreleased(:,t));
    ReleasedBSBLbs(k,t)=mean(BSBreleasedLbs(:,t));
    ReleasedBSBnumber(k,t)=mean(BSBreleased(:,t));
    ReleasedWLbs(k,t)=mean(WreleasedLbs(:,t));
    
    %Calculating total number of trips, aggregate catch and welfare:
    AggregateKeepSFLinear(k,t)=probTripLinear(k,t)*SFLbs(k,t);
    AggregateKeepSFnumberLinear(k,t)=probTripLinear(k,t)*SFnumber(k,t);
    AggregateReleaseSFLinear(k,t)=probTripLinear(k,t)*ReleasedSFLbs(k,t);
    AggregateReleasedSFnumberLinear(k,t)=probTripLinear(k,t)*ReleasedSFnumber(k,t);

    AggregateKeepSLinear(k,t)=probTripLinear(k,t)*SLbs(k,t);
    AggregateKeepSnumberLinear(k,t)=probTripLinear(k,t)*Snumber(k,t);
    AggregateReleaseSLinear(k,t)=probTripLinear(k,t)*ReleasedSLbs(k,t);
    AggregateReleasedSnumberLinear(k,t)=probTripLinear(k,t)*ReleasedSnumber(k,t);
    
    AggregateKeepBSBLinear(k,t)=probTripLinear(k,t)*BSBLbs(k,t);
    AggregateKeepBSBnumberLinear(k,t)=probTripLinear(k,t)*BSBnumber(k,t);
    AggregateReleaseBSBLinear(k,t)=probTripLinear(k,t)*ReleasedBSBLbs(k,t);
    AggregateReleasedBSBnumberLinear(k,t)=probTripLinear(k,t)*ReleasedBSBnumber(k,t);
   
    AggregateKeepWLinear(k,t)=probTripLinear(k,t)*WLbs(k,t);
    AggregateKeepWnumberLinear(k,t)=probTripLinear(k,t)*Wnumber(k,t);
    AggregateReleaseWLinear(k,t)=probTripLinear(k,t)*ReleasedWLbs(k,t);
    AggregateReleasedWnumberLinear(k,t)=probTripLinear(k,t)*ReleasedWnumber(k,t);
    
    TotalWelfareLinear(k,t)=log(exp(vLinear(k,t))+exp(opt_outLinear(k,t)))/(-parLinear(t,5)/1000); 
    NumberOftripsLinear(k,t)=probTripLinear(k,t);
    end
  end

k3=k;
  
%4. CARA Utility Model
    k=0;
  while (k<M) && (mean(sum(NumberOftrips_cara,1))<=TRIPS)
        k=k+1;
    %Random draw for trip non-catch characteristics
    IndvCharac=Y(randperm(end),:); 
    CostCharac=datasample(C,1); 
    CostCharacOpt_out=datasample(C,1);
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t);
    SEncounters=Sreshaffle(k,:,t); 
    BSBEncounters=BSBreshaffle(k,:,t);     
    WEncounters=Wreshaffle(k,:,t);
    
    for j=1:N   
     %Random draws for summer flounder
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
             
     %Random draws for scup
     if SEncounters(j)>0
     z=rand(1,SEncounters(j));
     Slegal=(z>=p_cutoff_S);
     Skept(j,t)=sum(Slegal,'all'); 
     Sreleased(j,t)=SEncounters(j)-Skept(j,t);       
     else
     Skept(j,t)=0;
     Sreleased(j,t)=0;
     end
     
     %Random draws for black sea bass
     if BSBEncounters(j)>0
     z=rand(1,BSBEncounters(j));
     BSBlegal=(z>=p_cutoff_BSB);
     BSBkept(j,t)=min(BSBbaglimit,sum(BSBlegal,'all')); %
     BSBreleased(j,t)=BSBEncounters(j)-BSBkept(j,t);       
     else
     BSBkept(j,t)=0;
     BSBreleased(j,t)=0;
     end
     
     %Random draws for weakfish
     if WEncounters(j)>0
     z=rand(1,WEncounters(j));
     Wlegal=(z>=p_cutoff_W);
     Wkept(j,t)=min(Wbaglimit,sum(Wlegal,'all')); 
     Wreleased(j,t)=WEncounters(j)-Wkept(j,t);       
     else
     Wkept(j,t)=0;
     Wreleased(j,t)=0;
     end   
     
     u_cara(t,j)=(1-exp(-par_cara(t,5)*(par_cara(t,2)*SFkept(j,t)/100+par_cara(t,1)*SFreleased(j,t)/100+par_cara(t,4)*(Skept(j,t)+BSBkept(j,t)+Wkept(j,t))/100+par_cara(t,3)*(Sreleased(j,t)+BSBreleased(j,t)+Wreleased(j,t))/100)))/par_cara(t,5);   
     end %End of the N random draws per individual
   
    vartheta_cara(t,:)=[par_cara(t,8),par_cara(t,9),par_cara(t,10),par_cara(t,11),par_cara(t,12),par_cara(t,13), par_cara(t,14)];
    
    %Calculating the individual utility for individual k of M
    v_cara(k,t)=mean(u_cara(t,:))+par_cara(t,6)*CostCharac(1); 
    opt_out_cara(k,t)=par_cara(t,7)+IndvCharac(1,:)*vartheta_cara(t,:)'+par_cara(t,6)*CostCharacOpt_out(1);
    
    %Calculating the probabilities of taking the trip
    probTripCARA(k,t)=exp(v_cara(k,t))/(exp(v_cara(k,t))+exp(opt_out_cara(k,t)));
    
    %Calculating catch per trip
    SFLbs(k,t)=mean(SFkeptLbs(:,t));
    SFnumber(k,t)=mean(SFkept(:,t));
    SLbs(k,t)=mean(SkeptLbs(:,t));
    Snumber(k,t)=mean(Skept(:,t));
    BSBLbs(k,t)=mean(BSBkeptLbs(:,t));
    BSBnumber(k,t)=mean(BSBkept(:,t));
    WLbs(k,t)=mean(WkeptLbs(:,t));
    Wnumber(k,t)=mean(Wkept(:,t));
    ReleasedSFLbs(k,t)=mean(SFreleasedLbs(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    ReleasedSLbs(k,t)=mean(SreleasedLbs(:,t));
    ReleasedSnumber(k,t)=mean(Sreleased(:,t));
    ReleasedBSBLbs(k,t)=mean(BSBreleasedLbs(:,t));
    ReleasedBSBnumber(k,t)=mean(BSBreleased(:,t));
    ReleasedWLbs(k,t)=mean(WreleasedLbs(:,t));
    
    %Calculating total number of trips, aggregate catch and welfare
    AggregateKeepSF_cara(k,t)=probTripCARA(k,t)*SFLbs(k,t);
    AggregateKeepSFnumber_cara(k,t)=probTripCARA(k,t)*SFnumber(k,t);
    AggregateReleaseSF_cara(k,t)=probTripCARA(k,t)*ReleasedSFLbs(k,t);
    AggregateReleasedSFnumber_cara(k,t)=probTripCARA(k,t)*ReleasedSFnumber(k,t);
    
    AggregateKeepS_cara(k,t)=probTripCARA(k,t)*SLbs(k,t);
    AggregateKeepSnumber_cara(k,t)=probTripCARA(k,t)*Snumber(k,t);
    AggregateReleaseS_cara(k,t)=probTripCARA(k,t)*ReleasedSLbs(k,t);
    AggregateReleasedSnumber_cara(k,t)=probTripCARA(k,t)*ReleasedSnumber(k,t);
    
    AggregateKeepBSB_cara(k,t)=probTripCARA(k,t)*BSBLbs(k);
    AggregateKeepBSBnumber_cara(k,t)=probTripCARA(k,t)*BSBnumber(k);
    AggregateReleaseBSB_cara(k,t)=probTripCARA(k,t)*ReleasedBSBLbs(k);
    AggregateReleasedBSBnumber_cara(k,t)=probTripCARA(k,t)*ReleasedBSBnumber(k);
        
    AggregateKeepW_cara(k,t)=probTripCARA(k,t)*WLbs(k);
    AggregateKeepWnumber_cara(k,t)=probTripCARA(k,t)*Wnumber(k);
    AggregateReleaseW_cara(k,t)=probTripCARA(k,t)*ReleasedWLbs(k);
    AggregateReleasedWnumber_cara(k,t)=probTripCARA(k,t)*ReleasedWnumber(k);
      
    TotalWelfare_cara(k,t)=log(exp(v_cara(k,t))+exp(opt_out_cara(k,t)))/(-par_cara(t,6)/1000); 
    NumberOftrips_cara(k,t)=probTripCARA(k,t);
    end
  end

k4=k;
  
varNames={'k_ambiguity','k_mixed','k_linear','k_CARA'};
writetable(table(k1,k2,k3,k4,'VariableNames',varNames),'k_NJ_2019.xls','Sheet',1);

toc;
save('CalibratedSeasonNJ2019.mat');
